#include "Collision1D_ChargeExchange.h"
#include "../Parallel/ParallelMPI.h"
#include "../Field/Field2D.h"


#include <cmath>
#include <iomanip>
#include <algorithm>
#include <ostream>
#include <sstream>

using namespace std;

Collision1D_ChargeExchange::Collision1D_ChargeExchange(PicParams* params_, CollisionParameter collision_param_, int collision_number_) :
    Collision1D(params_, collision_param_, collision_number_)
{
    cout<<"collision111"<<endl;
    read_cross_section();
    cout<<"collision222"<<endl;
    energy_threshold = crossSection[0][0];
    cout<<"collision333"<<endl;
}

Collision1D_ChargeExchange::~Collision1D_ChargeExchange()
{

}


void Collision1D_ChargeExchange::collide(PhysicalField* fields, vector<Species*>& vecSpecies, Diagnostic* diag, int itime)
{
    vector<int> index1, index2;
    vector<int> n1, n2;
    int n1_max, n2_max;
    vector<double> density1, density2;
    double density1_max, density2_max;
    vector<double> velocity_temp(3, 0.0), position_temp(1, 0.0);
    int collide_times = 0;
    int npairs_max;
    double npairs_double_max;
    unsigned int i1, i2, i3;
    double m1, m2, m3;

    double  sigma_cr, sigma_cr_max, ke1, ke_primary, ke_secondary,
            ran, P_collision;
    double  v_square, v_magnitude, v_magnitude_primary, v_magnitude_secondary;

    int i_space_global;
    double ke_radiative;


    if(collision_param.timestep_collision == 0 || itime % collision_param.timestep_collision != 0)
    {
        return;
    }


    Diagnostic1D *diag1D = static_cast<Diagnostic1D*>(diag);

    //electons
    auto s1 = static_cast<Species1D*>(vecSpecies[collision_param.species_number_list[0]]);

    //atom
    auto s2 = static_cast<Species1D*>(vecSpecies[collision_param.species_number_list[1]]);

    m1 = s1->species_param.mass;
    m2 = s2->species_param.mass;

    int n_space = params->dims_space[0];

    n1.resize(n_space);
    n2.resize(n_space);
    density2.resize(n_space);

    density2_max = 0.0;
    n1_max = 0;

    for(int i_space = 0; i_space < n_space; i_space++)
    {
        n1[i_space] = s1->particles_[i_space]->size();
        n2[i_space] = s2->particles_[i_space]->size();

        density2[i_space] = n2[i_space] * params->species_param[collision_param.species_number_list[1]].weight;

        if(density2[i_space] > density2_max) { density2_max = density2[i_space]; }
        if(n1[i_space] > n1_max) { n1_max = n1[i_space]; }
    }

    collide_times = 0;
    sigma_cr_max = max_cv(s1, s2);

    //npairs_double_max = n1_max * (1 - exp(-density2_max * sigma_cr_max * collision_param.timestep_collision * params->dt * collision_param.time_zoom_factor));

    for(int i_space = 0; i_space < n_space ; i_space++)
    {
        npairs_double_max = n1[i_space] * (1.0 - exp(-density2[i_space] * sigma_cr_max * collision_param.timestep_collision * params->dt));
        npairs_double_max *= collision_param.time_zoom_factor;

        //determine if zoom collision frequency
        if(collision_param.region_collision_zoom_start.size() > 0 && collision_param.region_collision_zoom_end.size() > 0)
        {
            if(params->local_min[0] + (i_space + 1) * params->cell_length[0] < collision_param.region_collision_zoom_start[0]
               || params->local_min[0] + i_space * params->cell_length[0] > collision_param.region_collision_zoom_end[0])
            {
                npairs_double_max /= collision_param.time_zoom_factor;
            }
        }

        if(n1[i_space] == 0 || n2[i_space] == 0)
        {
            continue;
        }

        auto cell_particles1 = s1->particles_[i_space];
        auto cell_particles2 = s2->particles_[i_space];

        index1.resize(n1[i_space]);
        random_shuffle(index1.begin(), index1.end());

        index2.resize(n2[i_space]);
        random_shuffle(index2.begin(), index2.end());

        npairs_double_[i_space] = npairs_double_max;

        npairs_[i_space] = npairs_double_[i_space];

        //See equations in http://dx.doi.org/10.1063/1.4742167
        npairs_rem_[i_space] += ( npairs_double_[i_space] - npairs_[i_space] );
        if(npairs_rem_[i_space] >= 1.0)
        {
            npairs_rem_[i_space] = npairs_rem_[i_space] - 1.0;
            npairs_[i_space]++;
        }

        if(npairs_[i_space] > n1[i_space])
        {
            cout<<"npairs is larger than the particle number in a cell!!!"<<endl;
            cout<<"npairs, n1 are: "<<npairs_[i_space]<<" "<<n1[i_space]<<endl;
            npairs_[i_space] = n1[i_space];
        }
        if(npairs_[i_space] > n2[i_space])
        {
            cout<<"npairs is larger than the particle number in a cell!!!"<<endl;
            cout<<"npairs, n2 are: "<<npairs_[i_space]<<" "<<n2[i_space]<<endl;
            npairs_[i_space] = n2[i_space];
        }


        for(int i = 0; i < npairs_[i_space]; i++)
        {
            //MESSAGE("nparis111"<<"  "<<i);
            i1 = index1[i];
            i2 = index2[i];
            auto& p1 = cell_particles1->data_[i1];
            auto& p2 = cell_particles2->data_[i2];

            v_square = (p1.velocity[0] - p2.velocity[0]) * (p1.velocity[0] - p2.velocity[0]) +
                       (p1.velocity[1] - p2.velocity[1]) * (p1.velocity[1] - p2.velocity[1]) +
                       (p1.velocity[2] - p2.velocity[2]) * (p1.velocity[2] - p2.velocity[2]) ;
            v_magnitude = sqrt(v_square);

            //kinetic energy of species1 (incident electrons)
            ke_primary = 0.5 * m1 * v_square / params->const_e;

            sigma_cr = v_magnitude * interplate_cross_section( ke_primary / params->const_e );
            P_collision = sigma_cr / sigma_cr_max;

            v_square = p1.velocity[0] * p1.velocity[0] + p1.velocity[1] * p1.velocity[1] + p1.velocity[2] * p1.velocity[2];
            v_magnitude = sqrt(v_square);
            ke_primary = 0.5 * m1 * v_square / params->const_e;

            //Generate a random number between 0 and 1
            double ran_p = (double)rand() / RAND_MAX;

            if(ran_p < P_collision)
            {
                calculate_scatter_velocity(p1, p2, m1, m2);

                velocity_temp[0] = p1.velocity[0];
                velocity_temp[1] = p1.velocity[1];
                velocity_temp[2] = p1.velocity[2];
                position_temp[0] = p1.position[0];

                p1.velocity[0] = p2.velocity[0];
                p1.velocity[1] = p2.velocity[1];
                p1.velocity[2] = p2.velocity[2];
                p1.position[0] = p2.position[0];

                p2.velocity[0] = velocity_temp[0];
                p2.velocity[1] = velocity_temp[1];
                p2.velocity[2] = velocity_temp[2];
                p2.position[0] = position_temp[0];

                collide_times++;

                v_square = p1.velocity[0] * p1.velocity[0] + p1.velocity[1] * p1.velocity[1] + p1.velocity[2] * p1.velocity[2];
                v_magnitude = sqrt(v_square);
                //kinetic energy of new species1 (ion)
                ke_secondary = 0.5 * m1 * v_square / params->const_e;
                ke_radiative = ke_primary - ke_secondary;
                i_space_global = params->local_min[0] / params->cell_length[0] + i_space;
                diag1D->radiative_energy_collision[collision_number][i_space_global] += ke_radiative;
            }
        }
    }
}



//elastic collision is assumed, same as Collisions1D_DSMC::ELASTIC
void Collision1D_ChargeExchange::calculate_scatter_velocity(Particle<1>& p1, Particle<1>& p2, double m1, double m2)
{
    double VRCP[3];		//VRCP(3) are the post-collision components of the relative velocity
    double VCCM[3];		//VCCM(3) are the components of the centre of mass velocity
    double RML;
    double RMM;
    double A,B,C,D;
    double OC,SC;
    double VR;

    double reduced_mass = m1 * m2 / (m1 + m2);
    RML = reduced_mass / m1;
    RMM = reduced_mass / m2;

    for (int i=0;i<3;i++)
    {
        VCCM[i] = RML * p1.velocity[i] + RMM * p2.velocity[i];
    }

    VR = sqrt( pow(p1.velocity[0] - p2.velocity[0], 2)
              +pow(p1.velocity[1] - p2.velocity[1], 2)
              +pow(p1.velocity[2] - p2.velocity[2], 2) );

    //use the VHS logic
    double RF = (double)random() / RAND_MAX;
    B = 2.0 * RF - 1.0;	//B is the cosine of a random elevation angle
    A = sqrt(1.0 - B * B);
    VRCP[0] = B * VR;

    RF = (double)random() / RAND_MAX;
    C = 2.0 * params->const_pi * RF; //C is a random azimuth angle
    VRCP[1] = A * cos(C) * VR;
    VRCP[2] = A * sin(C) * VR;

    for (int i = 0; i < 3; i++)
    {
        p1.velocity[i] = VCCM[i] + VRCP[i] * RMM;
        p2.velocity[i] = VCCM[i] - VRCP[i] * RML;
    }
}

